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I. Introduction 

This report contains user instructions, a listing and 
documentation for a microcomputer BASIC program that can be used 
to compute an estimate of the probability that a magnetometer 
based detection system such as the AN/ASQ-81 will detect a 
submarine magnetic dipole field during an encounter. 

The program generates detection probabilities based on two 
encounter models. In the first encounter model, the detection 
system uses a crosscorrelation detector. In the second encounter 
model, the detection system uses a square law detector. Relative 
to operationally realizable values, probabilities based on the 
first model represent upper bounds and those based on the second 
model represent lower bounds. For both encounter models, the 
signal is proportional to the square of the magnitude of the 
projection of the dipole field on the earth magnetic field and 
the magnetic noise does not change with changes in the position 
of the magnetometer. Also, for both encounter models, an 
encounter is a straight line encounter with constant vertical 
separation. 

The encounter models can be interpreted as models of a 
magnetic anomaly detection (MAD) system on an aircraft that is 
moving with constant course, speed and altitude in an encounter 
with a submarine moving with constant course, speed and depth. 

Or, they can be interpreted as models of a stationary magnetic 
anomaly detection system in an encounter with a submarine moving 
with constant course, speed and depth. 
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The program parameters include: encounter latitude and 

longitude, submarine induced magnetic moments, submarine 
permanent magnetic moments, submarine course, speed and altitude, 
magnetometer course, speed and depth, encounter lateral range 
(the horizontal range at the closest point of approach in a 
straight line encounter) and false alarm rate. (In the program, 
a false alarm is the event that the detection system classifies 
noise as a dipole signal.) 
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II. Program User Instructions 

As listed in Appendix 10, the program can be run under a 
BASIC language that is compatible with IBM PC BASIC. If the 
listing is used to enter the program through the keyboard, then 
the program should be saved with the name MAD. BAS or the value of 
N$ on line 40 should be changed to the file name under which it 
is saved. 

The program contains user instructions in the form of query 
and parameter limitation messages. As an example of the former, 
after starting the program under BASIC, the following message 
should appear: 

Magnetic Anomaly Detection (MAD) Lateral Range Function 

generate or print a program data file (g/p)? 

By entering p, data can be printed from a program data file 
that was generated by the program. By entering g, a program 
data file can be generated for a set of user specified 
conditions. With either response, a sequence of additional 
queries is displayed. These queries require either an indicated 
response or a parameter value as the input. If the initial 
response is g, the sequence includes queries whose responses 
determine whether or not an auxiliary data file will be generated 
that can be used for future input of magnetic, processing or 
kinematic data. In particular, the first query in the sequence 
gives the option of using a combined magnetic, processing and 
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kinematic data file. If the response indicates that it should be 
used and the file is available, the parameter values that remain 
to be entered in order to generate a program data file are the 
following: the false alarm rate, the magnetic noise, the maximum 
encounter lateral range and the lateral range step. The combined 
file should be used only if the effect of varying just one or 
more of these parameters is desired. If the response indicates 
that the file will not be used, gueries concerning magnetic, 
kinematic and processing parameter values are displayed. 

After all of the program parameter values have been entered, 
a query is displayed giving the option of generating the combined 
file. Then a query giving the option of generating a program 
data file, a query giving the option of printing the encounter 
parameter values and a query giving the option of printing 
lateral range function values are displayed. The lateral range 
function values are the encounter detection probabilities indexed 
by lateral range. The parameters maximum lateral range and 
lateral range step determine the index lateral ranges of the 
encounters for which probabilities of detection are computed. 

The program generates magnetic signal values that correspond 
to points in time during an encounter. Following the lateral 
range function query, a query is displayed that gives the option 
of printing the magnetic signal values for an encounter. If the 
option is exercised, the option is repeated. When the option is 
not exercised, a query is displayed that gives the option of 
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generating or printing a new program data file. If this option 
is not exercised, the program ends. 

Some suggested guides for determining parameter values can 
be found in Section IV of this report. 
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III. Encounter Model Limitations 



In the two encounter models, the signal is the unfiltered 
signal that is generated by a magnetic dipole that moves relative 
to a magnetometer. Depending on the input filtering, describing 
a ship magnetic anomaly as a dipole anomaly should not be a 
significant limitation for an encounter slant range at the 
closest point of approach (CPA) that is greater than one hull 
length. The detection decision is based on samples from a single 
time interval (window) that is centered on the CPA. The length 
of the sampling interval and the sampling rate are parameter 
values that are inputs to the program. In terms of signal-to- 
noise ratio, there is an optimum sampling interval length 
(integration time) and sample rate. Although dipole signal 
energy is not symmetrically distributed about the CPA time, for a 
given sampling interval, the difference between the signal energy 
for the optimum interval location and the CPA centered location 
should not be significant in most cases. 

The encounter model magnetometer noise samples are 
determined by a gaussian random process. They are values of 
independent identically distributed normal random variables. The 
standard deviation of these random variables is referred to as 
the magnetic noise and the variance is the noise in the sense of 
the signal-to-noise ratio. The degree of correspondence between 
this process and operational noise depends on the nature of the 
dominant operational noise sources and on the magnetometer input 
filtering (noise whitening) . 
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In the encounter models, the intervals are adjacent but not 
overlapping and a detection statistic corresponds to each sample 
interval and its value is determined by the sample values. If 
the value of the detection statistic for an interval equals or 
exceeds a threshold value, a detection is indicated. The 
threshold value is determined by the false alarm probability 
which in turn is determined by the false alarm rate and the 
sampling interval length. 

A detector that used a moving sample interval that was 
generated by replacing the oldest sample by the newest one would 
correspond more closely to the detector in an operational 
detection system. In a model of the detector, since the sample 
windows overlap, the detection statistics would represent a 
sequence of random variables that were correlated over an 
interval equal to the width of the sample interval. Because of 
this dependence, it seems unlikely that the results that would be 
obtained with an encounter model based on an overlapping interval 
detector would differ significantly from those obtained with an 
encounter model based on a nonoverlapping interval detector. 

For encounter lengths of the order of a few nautical miles 
or less, the straight line encounter condition should not be a 
significant limitation. In particular, this should be the case 
for a fixed magnetometer since, for a submarine (or surface ship) 
target, vertical separation and course changes should be less 
likely to occur. 
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Other models are available that can be used as the basis for 



computing an estimate of the probability that a magnetic anomaly 
detection system will detect a submarine during an encounter. 

For example, one that is described in Appendix 8 can be used to 
determine the slant ranges of straight line encounters for which 
the detection probability is equal to a specified value. The 
parameter values that are required to do this are an average 
submarine dipole moment, a detection system capability factor and 
a noise factor. Values for these parameters can be determined 
from operational data. However, the values are specific to 
averages over a particular set of encounter conditions. An 
advantage of the two encounter models relative to this model is 
their adaptability to different magnetic, processing and 
kinematic conditions. 
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IV. Parameter Values 



The magnetic parameter value queries are generally 
explanatory with regard to the value that should be entered. 

This is also true of the kinematic parameter value queries. 
However, there is some ambiguity with respect the processing 
parameter value queries and the noise parameter value query. To 
reduce this ambiguity, a brief discussion of the common 
characteristics of the two encounter model processing parameters 
and the noise parameter is given below. This is followed by some 
guidelines for choosing these parameter values. 

In both encounter models, a decision is made at the end of 
each sampling interval. The decision is either noise energy was 
present during the interval or noise energy plus signal energy 
was present during the interval. The sample intervals are 
adjacent, equal width, nonoverlapping time intervals. The number 
of samples that are input in a sample interval is determined by 
the sampling rate and the interval length. 

The program default choice for the sampling rate is 2 -MAXF 
where MAXF is a parameter that is labeled the maximum magnetic 
signal frequency. This sample rate is the Nyquist rate for an 
ideal low pass filter. However, the signal in a sample interval 
that is computed by the program represents an unfiltered dipole 
signal. This is a reasonable approximation if the signal energy 
that is associated with signal components greater than MAXF is 
relatively small. As discussed below, the noise energy in a 
sample interval should be considered to be proportional to MAXF 
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in order to be consistent with the encounter models. Ideally, a 
default choice for MAXF would make the ratio of the signal energy 
to the noise energy a maximum for the sampling interval of an 
encounter. The program default choice for MAXF is 2 -MAXVM/MINRO 
where MAXVM is a user estimated maximum encounter relative speed 
converted from knots to meters per second and MINRO is a user 
estimated maximum slant range at CPA in meters in terms of a just 
detectable target. (A precise definition of a just detectable 
target can be made in terms of a specified detection probability, 
false alarm probability and target dipole moment.) The default 
choice for MAXF is consistent with the observation in Reference 1 
that if an optimum value for MAXF is determined for a minimum 
dipole moment target, then no significant increase in MAXF is 
required in order to maintain a required detection probability if 
the encounter lateral range is decreased even though the signal 
energy spectrum is shifted to higher frequencies. 

Because of the detection statistics and the gaussian noise 
model, if there were no penalty for decision delay, a sample 
interval length for a signal should be chosen equal to the signal 
duration, since this would make the detection probability for an 
encounter a maximum. In the program, a default choice for the 
sample interval length (integration time) is 2 -MAXRO/MINVM where 
MINVM is a user estimated minimum encounter relative speed 
converted from knots to meters per second and MAXRO is a user 
estimated maximum slant range at CPA in meters for a detectable 
encounter (in terms of a specific detection probability and false 
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alarm probability) with a user specified maximum dipole moment 
target. This choice might be considered a balance between 
minimizing decision delay and maximizing detection probability. 

In the encounter models, the sample intervals are located so that 
the CPA time is at the center of a sample interval and this 
sample interval is the only one that contains signal energy. 

This characteristic is consistent with the default choice for the 
sample interval. For a given sample interval length, although in 
general the interval location is not the optimum one in terms of 
signal energy, it should be approximately so in most cases. 

The program noise parameter is SIG. It represents the 
standard deviation a associated with the magnetic noise process 
of the two encounter models. In terms of the ideal low pass 
filter implied by the encounter models, its square should be 
equal to MAXF-(SIGO) 2 where SIGO is the magnetic noise process 
constant power spectral density. The program does not enforce 
this relation. Therefore, in using the program, the implied 
relation between the two input parameters: magnetic noise and 
maximum magnetic signal frequency should be kept in mind. If an 
average value of the peak-to-peak magnetic noise for an encounter 
can be estimated, for example from a magnetometer trace, then the 
value for SIG should be chosen so that the estimate is 4 to 6 
times this value. 
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Appendix 1. The Detection Statistics 

In this appendix, y,, y 2 , • • • , y M are sequential values 

(voltages) representing the sample values in a sample interval. 
They are the input to a magnetometer's detector. With these M 
sample values, the detector computes the value of a detection 
statistic. This value is represented by x and the detection 
decision corresponding to the sample interval is determined by 
the decision rule: If x > x‘, then the input during the sample 
interval was noise plus signal, otherwise, the input was noise. 

For both encounter models, the detection probability and the 
false alarm probability are decreasing functions of x* and the 
relation is one-to-one in both cases. In the program, the false 
alarm probability p F is used to determine a unique value of 
the threshold x‘. This value is then used to determine a unique 
value of the detection probability p D . 

In the program, p F is found using the relation p F = R-<5t 
where R is the false alarm rate in false alarms per second 
and <5t is the sample interval length in seconds. This relation 
is based on the following argument: With no signal energy in a 
sample interval, y, = n, , y 2 = n 2 , • • • , y = n M where n,, n 2 , 

, n M are noise values (voltages) input to a magnetometer's 
detector that are values of independent normal (gaussian) random 
variables, each with mean zero and standard deviation a. 

Because of this, in the encounter models, values of x for 
different sample intervals are the values of independent random 
variables that determine two outcomes: x > x‘ or x < x*. 
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Therefore, p F is the same for each sample interval and, in 
terms of these outcomes, the sequence is a series of repeated 
independent Bernoulli trials and the expected number of trials 
between false alarms is l/p F . Since the time between trials is 
6t, the expected number of seconds between false alarms is 6t/p F 
and the expected number of false alarms per second is equal to 
p F /<5t. So, R = p F /<St. 

The determination of x‘ depends on the encounter model 
statistic. For both encounter models, when there is a signal, 

Yi ~ n i + s ,, y 2 — n 2 + y 2 , • , y M — n M + s M where s,, s 2 , • , 

s M are signal values (voltages) input to a magnetometer's 
detector. The models imply that the signal values s = c-H s 
where c is a constant whose value is determined by the 
characteristics of the encounter magnetometer and where the H s 
are dipole magnetic signal intensities. The models also imply 
that the noise values n are determined by a gaussian stochastic 
processes characterized by a standard deviation a and that 
n = c-H n + n' where the H N are magnetic noise intensities and 
the n' are magnetometer instrument noise values (voltages). In 
the program, the magnitude of c is 1. Since, for both models, 
p D depends only on the ratio of signal energy to noise energy 
for a sample interval, this is a satisfactory choice for the 
program if instrument noise process is assumed to be independent 
of the magnetic noise process and to be expressed in terms of an 
equivalent magnetic noise H = (1/c) n'. 
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For both encounter models, the signal (the average signal 
power) S = (1/M) -2 s K 2 and the noise (the expected value of the 
average noise power) N = a 2 so that the signal-to-noise ratio 
is (1/M) -2 s K 2 / cj 2 where the sum index K = 1, 2, ••• ,M. 

The Crosscorrelation Detector Statistic: The statistic for the 

first encounter model is a crosscorrelation detector statistic 
that is defined by the sum 

x = 2 y K -s K 

where the summation index K = 1, 2, • • • , M and the sum is over 

the values corresponding to a sample interval. For the first 
encounter model, the characteristics of both the noise and the 
signal are required in order to determine encounter detection 
probabilities. In particular, the signal values for an encounter 
are in the memory of the detector prior to the encounter. For 
the encounter conditions and a specified false alarm probability, 
the statistic is optimum in the sense that the encounter 
detection probability for this statistic is at least equal to 
that for any other statistic. Because of these considerations, 
encounter probabilities based on the crosscorrelation statistic 
can be considered to represent upper bounds on detection 
performance against dipole targets for magnetometers of the type 
described by the models. 

For a sample interval without signal energy, x is the value 
of a normal random variable with a mean ju x = 0 and a variance 
a x 2 = a 2 -2 s K 2 where a is the standard deviation associated with 
the noise process and the sum index K = 1, 2, • • • , M and the 



14 



sum is over the values corresponding to the sample interval. 

This implies that 

p F = 1 - P(x‘/ct x ) 

where P(z) is the standard normal cumulative distribution 
function. This relation is the basis for determining the 
threshold value x‘. 

For the sample interval with signal energy, x is the value 
of a normal random variable with a mean = E s K 2 where the sum 
index K = 1, 2, • • • , M and the sum is over the values 

corresponding to the sample interval. This implies that 

p D = 1 - P(v* - d 1/2 ) 

where v* = x'/a x and d = E s K 2 /a 2 = M-S/N. This relation is the 
basis for determining encounter detection probabilities for the 
first encounter model. The relation implies that for a specified 
false alarm probability p F the detection probability p D is an 
increasing function of the signal to-noise ratio S/N. 

The Square Law Detector Statistic: The statistic for the second 

encounter model is a square law (energy) detector statistic that 
is defined by the sum 

x = E y K 2 

where the sum index K = 1, 2, • • • , M and the sum is over the 

values corresponding to the sample interval. For the second 
encounter model, only the characteristics of the noise are 
required to determine encounter detection probabilities. 
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For a sample interval without signal energy, x/a 2 is the 
value of a chi-square random variable with M degrees of 
freedom. This implies that 

P F = 1 - P(x'/a 2 |M) 

where P(x*/a 2 |M) is the chi-square cumulative distribution 
function for a chi-square random variable with M degrees of 
freedom and where o is the standard deviation associated with 
the noise process. This relation is the basis for determining 
the threshold value x* . 

For the sample interval with signal energy, x/cr 2 is the 
value of a noncentral chi-square random variable with M degrees 
of freedom and noncentral parameter E s k 2 /ct 2 where the sum index 
K = 1, 2, • • • , M and the sum is over the values corresponding 

to the sample interval. This implies that 

p D = 1 - P ( x*/ a 2 1 M , E s K 2 / a 2 ) 

where P(x‘/a 2 |M, e s K 2 /a 2 ) is the cumulative distribution 
function for a noncentral chi-square random variable with M 
degrees of freedom and noncentral parameter E s K 2 /a 2 = M-S/N. 

This relation is the basis for determining encounter detection 
probabilities for the second encounter model. The relation 
implies that for a specified false alarm probability p F , the 
detection probability p D is an increasing function of the 
signal-to-noise ratio S/N. This is made more evident by the 
following relation: 

P(x‘/cr 2 |M, s s K 2 /a 2 ) = E{ (a J /J ! ) ■ exp (-a) • P[x*/a z | (M + 2 • J) ] } 
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where the parameter a = (1/2) -2 s */o 2 = (M/2) • (S/N) , the sum 



index J = 0, 1, 2, • • • and the sum index 

Note, since P(x‘/a 2 |M) < P[x‘/a 2 | (M + 2 • J) ] 

, the relation p D > p F is satisfied. 



K = 1, 2, • • • , M. 

for J = 0, 1, 2, 
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Appendix 2. Program Probability Calculations 

The program evaluates the cumulative and inverse cumulative 
distribution functions using approximations described in 
Reference 2. These approximations are listed below. 

The Standard Normal Cumulative Distribution Function 
Approximation : 

P(z) = 1 - s • t • (b, + t-(b 2 + t-(b 3 + t-(b 4 + t -b 5 ) ) ) ) 
where s = (1/ (2 - tt ) 1/2 ) -exp(-z 2 /2) and t = 1/(1 + b 0 -z). And 
where 

b 0 = .2316419, b, = .319381530, b 2 = -.356563782, 
b 4 = -1.821255978, b 5 = 1.330274429, 
and z > 0. For z < 0, P(z) = 1 - P(|z|). 

The Inverse Standard Normal Cumulative Distribution Function 
Approximation : 

Z (P) = t - (c 0 + t • (c, + t • c 2 ) ) / ( 1 + t-(d, + t • (d 2 + t-d 3 ))) 
where t = (In ( 1/Q 2 ) ) 1/2 and Q = 1 - P. And where 
c 0 = 2.515517, c, = .802853, c 2 = .010328, 
d, = 1.432788, d 2 = .189269, d 3 = .001308, 
and .5 < P < 1. For 0 < P < .5, z(P) = -z(l - P) . 

The Inverse Chi-Square Cumulative Distribution Function 
Approximation : 

v(P|M) = M • [ 1 - 2/(9-M) + z • (2/(9 -M) )' /2 ] 3 
where P(z) = P(v|M). In the program, the inverse standard 
normal cumulative distribution function approximation is used to 
determine z. 
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The Noncentral Chi-Square Cumulative Distribution Function 
Approximation : 

P(w|M,S s */o 2 ) = P(z) 

where z = [ 2 -w/ ( 1 + b) ] 1/2 - [ 2 • a/ ( 1 + b) - 1 ] 1/2 with 
a = M + S s k 2 /ct 2 , b = (2 s K 2 /a 2 )/(M + 2 s K 2 /a 2 ) and the sum index 

K=l,2, • • • , M. In the program, the standard normal 

cumulative distribution function approximation is used to 
determine P(z) . 



19 



Appendix 3. The Magnetic Signal 

The encounter models are defined by the following 
conditions: A submarine magnetic anomaly field is a magnetic 

dipole field that is superimposed on an earth magnetic field that 
is constant over an encounter region. A magnetometer magnetic 
signal value is the magnitude of the projection of a dipole 
magnetic field on the earth magnetic field at the location of the 
magnetometer. The basis for determining magnetic signal values 
in the program is an expression that can be developed as 
follows: In the right handed rectangular coordinate system that 

is shown in Figure 1, a magnetic dipole is at the origin, the xy- 
plane is the horizontal plane at a representative point in an 
encounter region, the positive y-axis is in the direction of 
magnetic north, the positive x-axis is in the direction of 
magnetic east and the positive z-axis is positive upward. In 
this rectangular coordinate system, the constant earth magnetic 
field can be expressed by H E = H E -(j-cos $ + k-sin $) where $ 
is a magnetic dip angle and H E is a magnetic field magnitude 
that characterizes the earth field in an encounter region. 

In a spherical coordinate system with the origin at the 
magnetic dipole and the polar axis in the direction of the dipole 
moment, H d = (c-p/r 3 ) -(2-r-cos 8 + 0-sin 8) is the magnetic 
field of the dipole at a point whose spherical coordinates are 
(r,r,8). In this expression p is the magnitude of the dipole 
moment and c is a constant whose value is determined by the 
choice of units. The magnetic signal for a magnetometer that is 
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Figure 1. A unit vector in the direction of the earth field and a 
unit vector in the direction of the dipole field of a 
dipole at the origin are shown at a point that is a 
distance r from the dipole. The unit vector h is 
in the direction of the horizontal component of the 
dipole moment. 
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described by the encounter models is H s = Hg-HyHg when the 
magnetometer is at the point (r,r,0) and H s can be expressed in 
terms of the rectangular coordinate system as follows: First, 

let r 0 be a unit vector in the direction of the spherical 
coordinate system polar axis. Then the magnetic dipole field 
H d = (c-p/r 3 ) • (3 -r -cos 6 - r 0 ) 

since r = (i-x + j-y + k-z)/r where r = (x 2 + y 2 + z 2 ) 1/2 , 

0 = [ (r 0 x r) x r]/sin 6 = [(r*r 0 ) -r - (r-r) -r 0 ]/sin e, and the 
dot product r*r 0 = cos 6. The unit vector r 0 can be expressed 
in the rectangular coordinates by noting that r 0 = p/p and then 
expressing p in rectangular coordinates. To do this, let 12 
be the depression angle of p from the xy-plane (the horizontal 
plane) with 12 positive downward and let a be the direction of 
p relative to magnetic north. Then the magnetic dipole moment 
p = p-(h-cos 12 - k-sin 12) in terms of the unit vector k and 
the unit vector h = i-sin a + j-cos 12 which has the direction 
of the horizontal component of p. With these results, the unit 
vector r 0 = i-(cos 12 -sin a) + j • (cos 12 -cos a) - k*sin 12 and 
H s = (c -p/r 3 ) • [ (3/r) -cos 6 • (y -cos $ - z -sin §) 

- (cos $-cos 12 -cos a + sin $-sin 12)] 
where cos 0 = (1/r) ■ (X'cos 12 -sin a + y-cos 12 -cos a - z-sin 12) 
since cos 0 = r*r 0 . As can be seen from this expression, for a 
constant dipole moment magnitude and direction and a constant 
earth field magnitude and direction, the magnetic signal is only 
a function of the rectangular coordinates of the location of the 
magnetometer relative to the dipole. In the encounter models, 
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both of these conditions are satisfied. However, by allowing p, 
n, a, H s and § to vary, the expression for H s is applicable 
to more general encounter models. 
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Appendix 4. The Anderson Formulation 

In the encounter models, the magnetic signal H s at a 
sample point in a straight line encounter can be represented in a 
form described by Anderson in Reference 4. For convenience, the 
Anderson formulation is used in the program to determine values 
for H s . It can be developed as follows: In Figure 2, the 

primed rectangular coordinate system is superimposed on the 
rectangular coordinate system of Figure 1 so that the origin is 
coincident with the origin of that system. A magnetometer is in 
a straight line encounter with a magnetic dipole that is located 
at the origin of the combined system. The combined system moves 
with the magnetic dipole with the x'-axis oriented so that it is 
parallel to and in the direction <p of magnetometer's relative 
to the magnetic dipole and the z'-axis oriented so that it is 
directed toward and passes through the CPA on that track. Let 
1, m, and n be the direction cosines of the dipole moment p 
and 1,, m, and n, be the direction cosines of the earth 
magnetic field H E in the primed system. Then the unit vector 
r 0 = i' -1 + j ' -n + k' -m. And, for points on the relative track, 
x' = s', y' =0 and z' = R where R is the slant range of the 
dipole at CPA and where s' is the algebraic distance of the 
magnetometer from CPA on the relative track. (It is negative for 
points before CPA and positive for points after CPA.) This 
implies that magnetometer position vector in the moving 
coordinate system is r = i'-(s'/r) + k'*(R/r). 
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Figure 2. The dipole is at the origin in the unprimed and primed 
coordinate systems. In the primed coordinate system, 
the CPA is at (0, 0, R) and, for a time t during 
the encounter, the magnetometer is at [s'(t), 0, R] . 
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